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Abstract: Nekrasov's integral equation describing water waves of permanent 
form, determines the angle <fi (s) that the wave surface makes with the horizontal. 
The independent variable s is a suitably scaled velocity potential, evaluated at 
the free surface, with the origin corresponding to the crest of the wave. For all 
waves, except for amplitudes near the maximum, <p (s) satisfies the inequality 
10001 <tt/6. 

It has been shown numerically and analytically, that as the wave amplitude 
approaches its maximum, the maximum of \4>{s) | can exceed 7r/6 by about 
1% near the crest. Numerical evidence suggested that this occurs in a small 
boundary layer near the crest where \4>{s)\ rises rapidly from |</>(0) | = and 
oscillates about 7r/6, the number of oscillations increasing as the maximum 
amplitude is approached. 

McLeod derived, from Nekrasov's equation, the following integral equation 



oo 

If sin <f>(t) s-t 

cp(s) = — / ; log 

3W 1 + J* S in 0( T )dT 



s + t 



dt 



for </> (s) in the boundary layer, whose width tends to zero as the maximum wave 
is approached. He also conjectured that the asymptotic form of <fi (s) as s — > oo 
satisfies 

<j>(s) = ^ {l + As" 1 sin (/31og s + c) + o(s~ 1 )}, 
6 

where A, [3 and c are constants with (3 « • 71 satisfying the equation 

V3f3 tanh = 1. 

We solve McLeod's boundary layer equation numerically and verify the above 
asymptotic form. 
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1 Introduction 



This paper considers the numerical solution of the equation 



sin </> (t) 



SttJ 1 + f* sin (p (r) dr 



log 



s + 1 



s - t 



dt 



— /fc(M)W> (t)-iP(s)}dt, 



(1.1a) 
(1.1b) 



where 



tp (t) = log ( 1 + / sin (t) dr I and k (t, s) 



2s 



(1.2) 



This equation was derived by McLeod [1] to describe the boundary layer behav- 
ior of the solution , for large fi, near the origin of the equation 



sin <^ (t) 



3nJ m -i + f Q sin^ (r)dT 



log 



F(* + i) 



F(a-t) 



(1.3) 



where F (t) = sn (Kt/n) and sn denotes the Jacobian elliptic function with 
quarter periods K and iK' . Equation (1.3) was first formulated by Nckrasov [2] 
to describe waves of constant periodic form moving with constant speed on the 
surface of a non-viscous fluid that is either of infinite depth or on a horizontal 
bottom, when the flow is taken to be irrotational. The wave is assumed to be 
symmetric about its crest and the equation is derived by conformally mapping 
the the region of the flow under one wavelength onto the unit disc cut along 
the negative real axis. The generic point on the circumference of the disc is e ts , 
with — 7r < s < 7r, and s — corresponds to the crest. As the circumference 
is described in a clockwise direction from — tt to tt the horizontal coordinate 
decreases by one wavelength. Then the function (p^ is the angle that the wave 
surface makes with the horizontal. With this choice of coordinate M (s) is 
periodic with period 27r. For more details, see Nekrasov [2], [3] and [4] or Milne- 
Thompson [5]. The wave is assumed to be symmetric about its crest. Thus 
4*^ (s) is an odd 2ir periodic function of s with (p^ (0) = 0. The solution is 
unique provided the additional assumption, that the wave has only one peak 
and one trough per period, is made. This is 



(p^ (s) > 0, s G (0, tt) with </v (0) = 0„ (tt) = 0. 



(1.4) 



The constants K and iK' , the quarter periods of sn, are related to the depth h 
and wavelength, A, by the relation 



K'/K = h/X. 



(1.5) 
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As h — > oo we have K ^tt (K' — » oo) and F (i) — ► sin |f so that (1.3) is 
also applicable for infinite depth. Equation (1.1) is derived by writing s — s/i 
and writing </> M (sfi) — <p (s) and letting \x — > oo with s fixed. Then <j) (s) satisfies 
(1.1). The boundary layer behavior of the solution of (1.3) was established 
numerically by Chandler and Graham [6], who were able to obtain a solution 
with a maximum value of (p^ (s) — 30 • 3787 . . .° and to detect a small number 
of oscillations about (p^ = 30° for \i = 10 18 . 

The numerical difficulty posed by the boundary layer behavior of the solu- 
tions of (1.3) for large \x is over come, by Chandler and Graham [6], by using a 
non uniform mesh for the discretisation of (1.3). This consists of three regions: 
one to cope with the rapid variation of <j)^ (s) in the boundary layer, whose 
thickness is of order near the origin; a second to deal with the slower vari- 
ation away from the origin and a third for the transitional layer in between. 
For further references on the analytical properties of the solutions of (1.3) and 
related numerical results, see Chandler and Graham [6] and McLeod [1]. 

The purpose of this paper is to solve (1.1) numerically and show that the 
solution 4>{s) oscillates about <p(s) = tt/6 and obeys the formal asymptotic 
result of McLeod [1] that can be written in the form 

{°° ^4 1 
1 + sin (n/3 log s + c„) > as s -> oo, (1.6) 

n=0 S J 

where A n and C n are constants and (3 = • 71 . . . is the root of 

V3/3tanh Q 71 "/ 3 ) = l - ( L7 ) 

Equation (1.1) represents the solution in the boundary layer and can thus be 
solved with a uniform mesh size. However (1.1) has an additional complication 
compared with (1.3) in that the range of integration is infinite and the decay of 
the solution to its asymptotic limit is algebraic. This fact means that we require 
careful consideration in order to obtain an accurate numerical representation of 
the integral in (1.1). 
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2 The Numerical Method 



Following Chandler and Graham [6] we solve the integral equation in the form 
(1.1b). This formulation is better, for numerical purposes, because the integra- 
tion by parts that is used to convert (1.1a) to (1.1b), removes the logarithmic 
singularity, at t = s, which occurs in the kernel of (1.1a). Although the corre- 
sponding kernel of (1.1b) has a pole, the singularity of the integrand is removable 
since the multiple if) (t) — ip (s), has a simple zero at t = s. 
Thus we write 

<t>(s) = ^J o °°K(t,s)dt, (2.1) 

where 

2s (V; (t)-y>0Q) 

K(t,s) = g2 _ t2 t^s (2.2a) 

, / / n sm(b (t) .„ „, , 

= W = t-t* 1 w (2 ' 2b) 

1 + J sm (r) ctr 

the value in (2.2b) being the limit of the right hand side of (2.2a) as \t — s\ — > 0. 

We aim to set up a numerical approximation to the integral in terms of 
a discrete number of values <f>(si), where Sj = ih, < i < 2N, with ./V an 
integer, for suitable choices of h and N and a continuous set of values <fi (s) for 
s ^ 2Nh. Any values of <f> (s) for s < required by the numerical approximation 
are determined by the fact that (f> (s) is an odd function of s. The numerical 
representation of the integral requires two approaches. The first is a finite 
difference formulation of the integral over a predetermined finite range using 
the discrete values of <f> and the second is an estimation of the remainder using 
an appropriate asymptotic estimate of the values of 4> ( s ) f° r s ^ 2Nh. The 
details of the asymptotic form of <f> (s) as s — > oo that is used will be discussed 
later. 

So we choose an appropriate end point 2T where T is given by T = Nh and 
we can approximate the integral I\ (s, 4>) — J Q 2T K (t, s) dt using Simpson's Rule, 
since the integrand is analytic. The choice of the end point 2T is some what 
arbitrary. Eventually, see below, we will want to consider 1\ (s, 4>) for values 
of s < T. We choose an end point mT, with m = 2 in this case, so that the 
singularity of k (t, s) at t — s is far from the end point. The reason for this 
is that the remainder integral, again see below, requires a different evaluation 
and it is advantageous to make sure that the singularity of k(s,t) is not close 
to the range of t in the remainder integral. This will become clearer when the 
evaluation of the remainder integral is discussed later. 

Assuming that for large s, <f> (s) is known in the form of an asymptotic ex- 
pansion then truncation of this series, expansion of the integrand and a term 
by term integration of the integrand will give a suitable analytical estimate 
EI 2 (s, 4>) for the integral I 2 (s, (f>) = K (t, s) dt. Then we define the numcr- 
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ical representation of the integral in (2.1) as 

NI {s, <j>) = NI\ {s, <j>) + EI 2 (s, 4>) . (2.3) 

An alternative approach, assuming that the asymptotic form of 4>(s), s > T, 
has been chosen, is to transform the infinite range of the remainder integral into 
a finite range, which can then be approximated numerically. For this purpose 
it is more convenient to revert to the integral in the form (1.1a) so we write 

00 

h (s) = log (f^) (V> (2T) - $ (s)) + Jk 3 (s, t) dt, (2.4) 



2T 



where 



M«)= . ( *> log(^) (2-5) 



1 + J Q sin <p (t) dr 



t - s 



If 4>(t) — > 7r/6 + 0(i 1 ) and f^°(<f>(t) — ir/Q)dt is bounded, it is easily estab- 
lished that &3 (s,i) = 2st~ 2 + o(t~ 2 ) as t — > 00. Thus the integral of fc 3 , in 
(2.4) is convergent at infinity and the substitution i = 2T/u transforms it to 
J Q fc4 (s, u) du with &4 (s, 0) = s/T. This integral can now be approximated using 
Simpson's rule with a suitably chosen step length. This approximation can be 
used instead of EI2 (s,<p) in (2.3). 

Simpson's rule gives an approximation which is of order h , but this rule 
requires an interval which consists of an even number of step lengths. However 
the integrand contains the function ip (t) which involves the determination of 
J sin 4> (t) dr at values t = U = ih. To obtain a numerical approximation to this 
which is the same order as Simpson's rule for this integral we use an appropriate 
modified trapisoidal rule. 

We now wish to solve the approximation 

<M S ) = -^NI(s,4>)- (2.6) 

To do this we define an approximation (j>N (si) to the solution <j> (s) at the discrete 
values Si — ih, < i < N. Using the same asymptotic form at the solution as 
that used to define <j> (s) for s ^ 2Nh we define the remaining discrete values of 
(f)M {si), N + 1 < i < 2Nh, required for the evaluation of NIi (s) at the points 
s = Si, < i < N. 

Thus 4>n (si) satisfies the equations 

<Pn (Si) = ^-NI ( 8i , (j) N (sj)) ,0<i<N. (2.7) 

This gives, in a similar fashion to Chandler and Graham [6], a fully discrete 
non-linear system for the unknowns {4>n («i) , i — 0..N}. This system is solved 
by the iterative method 

<ffi ( Sl ) = NI (si,^™-^ ( 8j j) , i = 0..N, (2.8) 
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starting from a suitable initial approximation (jr N (sj). Chandler and Graham 
[6] were able to prove that, when the quadrature method used to approximate 
their integrals was the trapisoidal rule, convergence was guaranteed, although 
for computational purposes they opted for a more accurate scheme for computa- 
tional purposes. Their proof cannot be extended to the numerical approximation 
used here even if the quadrature method is the trapisoidal rule because of the 
infinite range of integration. However we find that, as in the cases looked at by 
Chandler and Graham [6], the convergence rule is very quick. 
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3 The necessity of rescaling 



We see from the definition of K (t, s) , (2.2a,b), and the fact that (j) (0) is zero, 
that NI (0, (f>) = provided the initial guess ^ (0) = 0. Then (2.7) gives 
^TV (0) = f° r au rn > 0. Thus effectively we can work with the N variables 
{</>jv (si) , i = 1..N} and corresponding N equations from (2.7). One of the aims 
is to verify the asymptotic result (1.6). Initially we do not assume this and 
report here that for a variety of sensible choices of the asymptotic form of <f> (s) 
we get rapid convergence to the solution of (2.8). Provided T is sufficiently 
large we can then numerically verify that (1.6) is the correct asymptotic result, 
using the computed values of <ft (s) for s < T. Having verified this numerically 
to get the best accuracy we use (1.6) and find that as well as providing a more 
accurate numerical solution the convergence rate is also improved. The larger 
T is, the less necessary it is to have a large number of terms from (1.6) and in 
practice we use 

cj>(s) = ^(l + ^ sinologs + c)j , s>T. (3.1) 

Table 1 shows the comparison of the location and the values of <j> (s) at successive 
maximum and minimum values of <j> and the comparison between this method at 
that of Chandler and Graham [6] . Before discussing this comparison we use the 
values of s at the successive turning points to illustrate the need for rescaling 
the variable s. It will become clear that the computations done to obtain table 
1 could not be achieved by the method outlined in paragraph 1. We see that 
the s coordinate of each successive turning point increases by a factor of about 
81, which is approximately the value of . This is compatible with the set 
of turning points obtained from (3.1). The last turning point in < s < T is 
located at s = 2 x 10 11 . Typically we used h = 1/20 as a sensible choice of h 
compatible with having a large enough T to capture the asymptotic behavior of 
the solutions. However with this choice of h it is not feasible to take T = 2 x 10 11 
as this would involve 4 x 10 12 grid points. Typically using the scheme outlined in 
paragraph 1 we chose T = 100 and this does not even get to the first minimum 
of 4> (s). However we learn from this initial attempt at a numerical solution that 
beyond s — 100, 6 \(f> (s) — n/6\ /n < 10~ 2 and varies very slowly. Thus for large 
s we do not need to take such a small step length. 

For the numerical scheme we have used, we require a constant steplength 
so we make a simple change of independent variable. We wish to make no 
effective change at the origin but an exponential change at infinity so we use 
the transformation s = e v — 1. Then with t = e z — 1 and 9 (y) — <f> (s (y j) , (1.1) 
becomes 

1 f°° 2 (e v — 1) 

We are then able to reduce the step length, h, and still take T = e VT — 1 to be 
large. Typically we take h = 1/100 and y T = 30 giving T = 1.0 x 10 13 . This 
requires 3000 unknowns <j> (Vi) where yi = ih, i = 1 .. 3000. 



1 + / Z sin 9(Q 
1 + sin 9(0 



dy, (3.3) 
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After the rescaling, the numerical scheme is essentially the same as that given 
in section 2 and is not repeated. However near y = dt, 6 \ 6 (y) — 7r/6| /it is now 
of order 1CP 13 so the form of 9 (y) effectively given by (3.1) will be accurate to 
1CT 26 , that is 0(T- 2 ). 
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4 The Numerical Results and Conclusions 



All the numerical results given here are those produced by the numerical scheme 
outlined in Section 2 and 3 using the rescaled problem. Table 1 shows the 
comparison of the successive maxima and minima of <f> (s) compared with those 
computed for the full problem by Chandler and Graham [6]. The position 
of these maxima and minima for the Chandler and Graham [6] computation, 
has been calculated by scaling their coordinate, s, by /i compatible with the 
boundary layer scaling used to derive (1.1) from (1.3). Thus s = sb-s = s CeG x 
ji. The number of decimal places given in table 1 for this numerical computation 
are as accurate as the numerical calculation will allow. There are three forms 
of error: the first comes from the order of the numerical approximation to the 
solution which is O (/i 4 ) which gives rise to errors of order 10 -8 ; the second is 
due to machine accuracy which gives rise to an error of about 10 -14 to 10~ 16 ; 
thirdly there is the error that arises when predicting the position and size of 
the maxima and minima of a function, from discrete data at given grid points, 
assuming that the data is accurate. The figures quoted in table 1 do not take 
into account the first of two of these sources of error. 

The comparison with the computations of Chandler and Graham [6] is very 
good. The value at the first maximum is the same to eight significant figures 
and the position the same to six significant figures. The calculation of the 
value at the maximum always being more accurate that its positions. The 
values at the first minimum are in similar agreement although Chandler and 
Graham [6] only quote the position to four significant figures and the value 
at the minimum is only 4 x 10~ 3 below 30° so relatively the numbers do not 
appear to be in such good agreement as the value at the first maximum. The first 
noticeable divergence of the two computations appears at the second minimum 
where the estimates of the positions differ by about 4% although the values at 
this minimum are in good agreement given that they are both of order 10~ 7 
below 30°. However the next maximum of Chandler and Graham [6] lies below 
30° and it is apparent that at this value of s the effects of the outer solution, 
that is the decrease from the maximum on a slower scale, are just beginning 
to show. Presumably at this value of /j, the oscillations in the Chandler and 
Graham [6] begin to cease at or around this value of s. 

We wish to show that the solution behaves like (1.6) for large s. So for 
comparison we write 6 (x) = <j> (s), where x = & logs so that we expect 

9 (x) <~ — < 1 + — sin7r (x — x ) + . . . 1 as x — > +oo (4.1) 
6 s 



or 



* (x) = I -6 (x) - 1 I s ~ Asimv (x - x ) + ■ ■ ■ , (4.2) 



7T 



Compared with the transformation (3.1) which has y — when s = we have 
x — > — oo as s — > 0. This makes (^6 (x) — l) s — > as x — > — oo and introduces 
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a minimum of the function \l/ (a;) before the first maximum. The values of 
x = Xi at the minima, maxima and the zeros of '5 (x) and the value of (x) 
at the turning points are shown in table 2. If (4.2) were to be exact then the 
difference Xi — Xi-\ — 1/2 = Axi would be zero and the magnitude of the value 
of $ (x) at the turning points would be constant and equal to A. Included in 
this table are the computed values of Axi . 

From the table we see that a good fit is obtained by choosing A and xq so 
that VE" (x) and (4.2) agree at the second maximum and fourth zero this gives 

A = 1 • 2364860386 ... and r = • 72422 .... (4.3) 

A plot of the asymptotic expression (4.2) with these values of A and Xq and the 
comparison with '3/ (x) is given in figure 1. The two graphs are indistinguishable 
from each other over a surprisingly large range of values of x, from before the 
first zero to beyond the sixth zero. The graphs start to diverge after this point. 
This is due to the fact that the exact solution of <fi (s) — 7r/6, or equivalently 
^(x)/s, is so small in this range that round off error starts to become important 
and eventually dominates the solution. This is more apparent in figures 2 and 3 
which plot the difference between 'J (x) and its asymptotic value. Figure 2 shows 
this difference multiplied by 100 in the range of values of x where the difference 
is less then one, while figure 3 shows 1000 times the difference. In both figures 
we see that the difference increases rapidly after x — 4. It is particularly visible 
in figure 3 that this rapid rise has two different components: a systematic rise 
due to truncation error of the numerical scheme, which is of order 10 8 and a 
random error on the scale of about 10 , due to machine accuracy. 

The last plot, figure 4, shows the difference between * (x) and its asymp- 
totic value multiplied by s. This clearly shows that the dominant feature is 
one of a periodic function of period 1, compatible with a term proportional to 
s~ 2 sin27r (x — xi) that appears in (1.6). 

To conclude we have presented a numerical scheme for the solution of (1.1), 
written in the form (3.3) which allows a sufficiently accurate numerical solution 
over a range < s ^ 10 13 , that we can verify the predicted asymptotic form 
(1.6). The numerical calculation is limited by the two factors, truncation error 
and machine accuracy. The numerical solutions can be made more accurate 
by a higher order integration scheme but the range of integration is limited 
because the difference between the solution and ir/6 becomes the same order of 
magnitude as the machine accuracy. 



10 



References 

1. J.B. McLeod, The Stokes and Krasovskii Conjectures for the wave of greatest 
height. Stud. App. Math. 98: 311-333 (1997) 



2. A.I. Nekrasov, Izv. Ivanovo- Vosnosonk. Politehn Inst. 3: 52-65 1921; 6:155- 
71 (1922) 

3. A.I. Nekrasov, Izv. Ivanovo- Vosnosonk. Politehn Inst. 6:155-71 (1922) 

4. A.I. Nekrasov, The exact theory of steady state waves on the surface of 
a heavy liquid. Technical Summary Report No 813. Mathematical Research 
center, University of Wisconsin, 1967 [D.V. Thampuran, translator:C. W. Cryer, 
editor] 

5. L.M. Milne-Thompson, Theoretical Hydrodynamics, Macmillan, London, 
1968. 

6. G.A. Chandler and I.G. Graham, The Computation of water waves modelled 
by Nekrasov's Equation. SI AM J. Numer. Anal. 30: 1041-1065 (1993). 



11 



Figure Captions 

Tabic 1. Positions of the turning points, s t and the corresponding values, 4>(st) 
and comparison with those obtained by Chandler and Graham. 
Table 2. The positions, Xi of the zeros and the turning points of s((j)(s) — n/6) 
as a function of x = f3logs and the corresponding values at the turning points. 
Axi is the difference Xi — Xi-\ — \ 

Figure 1. Comparison ^(x) = (69(x)/7r— l)s with Asin(7r(x — x )) as a function 
of x — f3 'logs /it. 

Figure 2. Difference between the solution and its Asymptotic form 100(^(x) — 
j4sin(7r(a; — .to))) as a function of x — /31ogs/7r. 

Figure 3. Difference between the solution and its Asymptotic form 10000(4' (x) — 
Asin(7r(x — x ))) as a function of x = /3\ogs/ir. 

Figure 4. Difference between the solution and its Asymptotic form 'i'i(x) = 
s((6Q(x)/tt — l)s — Asin(7r(x — xq))) as a function of x = /3\ogs/n. 
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Figure 1. Comparison #(x) = (6Q(x)/tt — l)s with Asin(7r(izi — x )) as a function of x = (3\og 
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Figure 2. Difference between the solution and its Asymptotic form 

100(^(a;) — ^4sin(7r(a; — x ))) as a function of x = f3\ogs/ir. 
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Figure 3. Difference between the solution and its Asymptotic form 

10000(^(0;) — Asin(-7r(a; — x ))) as a function of x = /3\ogs/w. 
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Figure 4. Difference between the solution and its Asymptotic form 

^>i(x) = s((6Q(x)/tt — l)s — Asin(7r(x — x )j) as a function of x = f3\ogs/n. 
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TABLE 1: Positions of the turning points, s t and the corresponding 
values, 4>(s t ) and comparison with those obtained by Chandler and 

Graham. 



■St This paper 


s t C. & G. 


60(Si)/7T - 1 


<f>{a t )° - 30° 


</>(s t )° -30° C. & G. 


5.706256X10 1 


5.7062493X10 1 


1.26234416 xl0~ 2 


3.787032480 xlO" 1 


3.787032466X10- 1 


4.683476245 xlO 3 


4.683xl0 3 


-1.5345108772 xlO" 4 


-4.6035326316xl0- 3 


-4.60353xl0~ 3 


3.80716.261 xlO 5 


3.807xl0 5 


1. 88776874 xlO" 6 


5.66330622 xlO" 5 


5.6631 xlO" 5 


3.09513xl0 7 


3.21xl0 7 


-2.322037xl0~ 8 


-6.966111xl0- 7 


-7.4218xl0~ 7 


2.51266xl0 9 


2.416xl0 8 


2.8545xl0~ 2 


8.5635xl0~ 9 


-3.6722xl0~ 7 


2.058xlO n 




-3.96xl0~ 12 


-1.188xl0~ 10 





TABLE 2: The positions, X{ of the zeros and the turning points of 
s((/)(s) — 7r/6) as a function of x = /3logs and the corresponding 
values at the turning points. Axi is the difference x% — Xi-\ — \ 
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